Onset of thermal convection in a horizontal layer of granular gas 
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The Navier-Stokes granular hydrodynamics is employed for determining the threshold of ther- 
mal convection in an infinite horizontal layer of granular gas. The dependence of the convection 
threshold, in terms of the inelasticity of particle collisions, on the Froude and Knudsen numbers is 
found. A simple necessary condition for convection is formulated in terms of the Schwarzschild's 
criterion, well-known in thermal convection of (compressible) classical fluids. The morphology of 
convection cells at the onset is determined. At large Froude numbers, the Froude number drops out 
of the problem. As the Froude number goes to zero, the convection instability turns into a recently 
discovered phase separation instability. 

PACS numbers: 45.70.Qj 



INTRODUCTION 

Fluidized granular media exhibit a plethora of fas- 
cinating pattern-formation phenomena that have been 
subjects of much recent interest In this work we 
address thermal (buoyancy-driven) granular convection 
1^, ^ ^ H, 1^. Being unrelated to the shear or time- 
dependence introduced by the system boundaries, it re- 
sembles the Rayleigh-Benard convection in classical fluid 
and its compressible modifications ||, |ll|, |l^ . In 
classical fluid convection requires an externally imposed 
negative temperature gradient, that is a temperature gra- 
dient in the direction opposite to gravity. In a vibroflu- 
idized granular medium a negative temperature gradient 
sets in spontaneously because of the energy loss by in- 
elastic collisions. Convection develops when the absolute 
value of the temperature gradient is large enough. In the 
simplest model of inelastic hard spheres that we will use 
it happens when the inelasticity coefficient q = (1 — r)/2 
exceeds a critical value depending on the rest of the pa- 
rameters of the system. Here r is the coefficient of normal 
restitution of particle collisions. 

Thermal granular convection was first observed in 
molecular dynamics (MD) simulations of a system of 
inelastically colliding disks in a two-dimensional (2D) 
square box ||]. The boundaries of the box ||] did not 
introduce any shear or time-dependence: the system was 
driven by a stress-free thermalizing base. The top wall 
was perfectly elastic, while the lateral boundaries were 
either elastic or periodic. Experiment with a highly flu- 
idized three-dimensional granular flow Q gives strong ev- 
idence for thermal convection, though energy loss at the 
side walls introduces complications ^. A clear iden- 
tification of thermal convection in experiment requires a 
large aspect ratio in the horizontal direction, so that mul- 
tiple convection cells can be observed. MD simulations 
in 2D of a vibrofluidized granular system with a large 
aspect ratio indeed show multiple convection cells @]. 

This work deals with a theory of thermal granular 
convection in a system with a large aspect ratio. Re- 
cently, a continuum model of thermal granular convection 



has been formulated g in the framework of the Navier- 
Stokes granular hydrodynamics. In the dilute limit the 
Navier-Stokes hydrodynamics (or, more precisely, gasdy- 
namics) is systematically derivable from more fundamen- 
tal kinetic equations jl^, |lj]. Like any other hydrody- 
namic approach, the Navier-Stokes hydrodynamics de- 
mands small Knudsen numbers for its validity. In addi- 
tion, it has been shown that at moderate inelasticities q 
non-hydrodynamic effects (such as the lack of scale sepa- 
ration, the normal stress difference and non-Gaussianity 
in the particle velocity distribution) may become impor- 
tant [|l5j. Therefore, the Navier-Stokes granular hydro- 
dynamics is expected to be accurate quantitatively only 
for nearly elastic collisions, q 1. Though restrictive, 
the nearly elastic limit is conceptually important. Also, 
one can expect some of the results, obtained in this limit, 
to be still qualitatively valid for larger inelasticities, such 
as those encountered in experiment. 

In Ref . Q the full set of nonlinear hydrodynamic equa- 
tions for thermal granular convection was solved numer- 
ically in a 2D box with aspect ratio 1. It was observed, 
in qualitative agreement with MD simulations that 
the static state of the system (a steady state with a zero 
mean flow) gives way to convection via a supercritical 
bifurcation, the inelasticity q being the control parame- 
ter. The present work employs the same hydrodynamic 
formulation |^ to perform a systematic linear stability 
analysis of the static state. We determine the convec- 
tion threshold as a function of the scaled parameters of 
the problem and of the horizontal wave number of small 
perturbations. This analysis makes it possible to predict 
the convection threshold and determine the morphology 
of the convection cells in a system with an arbitrary as- 
pect ratio, including an infinite horizontal layer, a stan- 
dard setting for convection in classical fiuids |Q, We 
also formulate a simple necessary (but not sufficient) cri- 
terion for thermal granular convection in terms of the 
Schwarzschild's criterion, well-known in thermal convec- 
tion of (compressible) classical fluids . Finally, we take 
the limit of a zero gravity and establish the connection 
between thermal convection and a recently discovered 
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phase-separation instability 0, |T|, |T|, |l|, |2|. 

MODEL AND STATIC STATE 

Let a big number of identical smooth hard disks with 
diameter d and mass m move and inelastically collide 
inside an infinite two-dimensional horizontal layer with 
height H. The gravity acceleration g is in the negative 
y direction. The system is driven by a rapidly vibrat- 
ing base. We shall model it in a simplified way by pre- 
scribing a constant granular temperature Tq at y = 0. 
The top wall is assumed elastic. Hydrodynamics deals 
with coarse-grained fields: the number density of grains 
n(r,t), granular temperature T{r,t) and mean flow ve- 
locity v(r,t). In the dilute limit, the scaled governing 
equations are 

dn/dt + nV -v = 0, (1) 
nd\/dt^y -P - Fney, (2) 

ndT/dt + nTV-v = KV ■ {T^/^VT) - KRn^ T^'^ . (3) 

Here d/dt — d/dt v • V is the total derivative, 
P = -nTl + (l/2)Xri/2f) is i-jjg g^j.ggg tensor, D = 
(1/2) (Vt; + (Vtj)-'") is the rate of deformation tensor, 
D = D — (1/2) tr (D ) I is the deviatoric part of D, and 
I is the identity tensor. In the dilute limit, the bulk 
viscosity of the gas is negligible compared to the shear 
viscosity , so only the shear viscosity is taken into ac- 
count. In addition, we have neglected the small viscous 
heating term in the heat balance equation (^. The in- 
elastic contribution to the heat flux is proportional 
to q and can be safely neglected at small q. In the 2D 
geometry, the three scaled parameters entering Eqs. (|^) 
and (H) are the Froude number F = mgH/To, the Knud- 
sen number K = 2t:~^^'^ {dH{n))~^ and the relative heat 
loss parameter R — 8qK~'^. Furthermore, (n) is the to- 
tal number of particles per unit length in the horizontal 
direction, divided by the layer height H. It will be con- 
venient to use the relative heat loss number R instead of 
q. In Eqs. (||)-(||), the distance is measured in the units 

1/2 

of H, the time in units of H/Tq , the density in units of 
(n), the temperature in units of Tq, and the velocity in 

1/2 

units of Tq . The Navier-Stokes hydrodynamic model 
is expected to be valid when the mean free path of the 
particles is much smaller than any length scale (and the 
mean collision time is much smaller than any time scale) 
described hydrodynamically. This implies, in particular, 
that the Knudsen number should be small: K ^ 1. 

The boundary conditions for the temperature are 
T{x, y = 0, t) = 1 at the base and a zero normal com- 
ponent of the heat flux at the upper wall: dT/dy{x, y = 



l,t) = 0. For velocities we demand zero normal compo- 
nents and slip (no stress) conditions at the boundaries. 
The total number of particles in the system is conserved: 

1™ ^ / '^^ / d.yn{x,yA) = 1, (4) 

L->oo ZIj J_i^ Jq 

where we have introduced the horizontal dimension 2L. 
The hydrodynamic problem is characterized by the three 
scaled numbers F, K and R. 

The simplest steady state of the system is the static 
state: no mean flow. At a nonzero gravity the density 
and temperature of the static state depend only on y, 
and are described by the equations 

{nsTs)' + Fns = and (T^/^T^)' - RnlT^/^ ^0, (5) 

where the primes denote the y-derivatives. The boundary 
conditions are Ts(0) — 1 and (1) ~ 0, the normahza- 
tion condition is dy Us (y) = 1 . The static state (see 
Fig. 1^) is characterized by two scaled numbers: F and 
_R, and can be found analytically by transforming from 
the Eulerian coordinate y to the Lagrangian mass coor- 
dinate ii{y) = ns{y')dy' ||]. At large enough R, a 
density inversion develops: a denser (heavier) gas is lo- 
cated on the top of an under-dense (lighter) gas. This 
is clearly a destabilizing effect that drives thermal con- 
vection. However, this effect is neither sufficient, nor 
necessary for convection, see below. The actual (nec- 
essary and sufficient) criterion should take into account 
the heat conduction and viscosity that scale like K in the 
governing equations and have a stabilizing role. A small 
sinusoidal perturbation in the horizontal direction is un- 
stable with respect to convection if the relative heat loss 
parameter R exceeds a critical value that depends on F, 
K and the horizontal wave number. 
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FIG. 1: Static temperature (solid line) and density (dashed 
line) profiles for _F = 0.1 and R — 0.5. 



THE LINEAR STABILITY ANALYSIS 

The linear stability analysis involves linearization of 
Eqs. (|l|)-(^) around the static solutions ns{y) and Ts{y). 
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The linearized equations are 



1/2 

We have denoted for brevity ao = Tg (y) , ai 



dh dvx d , N „ 



dt 



-FflBy, (7) 



dT 



where n, T and v denote small perturbations. Exploiting 
the translational symmetry of the static state in the hor- 
izontal direction, one can consider a single Fourier mode 



h{x,y,t)=e '^^ N{y) cos k^x , 



T{x,y,t) 



-It p. 



Q{y) cos k^x, 



Vx{x,y,t) = e ^*u(y)sinfca:X, 



Vy{x,y,t) 



v{y) cos kxX . 



(9) 



where 7 is the scaled growth/decay rate, and k^ is the 
scaled horizontal wave number. Substituting Eq. (|^) 
into Eqs. (^)-(^) and eliminating N, we obtain three 
homogeneous ordinary differential equations that can be 
written as a single equation for the eigenvector U(j/) = 
[Q{y), u{y), v{y)], corresponding to the eigenvalue 7: 



AU" + BU' + CU = 0, 



where 



Kaa 
Kao/4: 
Kao/4 -n,T,/7 

2Kai -2KRnlTs''^ h - ^sTs 

B = I Kai/4 k,n,Tjj 

-Us -k^UsTs/j K ai/A - Tsu'Jj 

while the elements of matrix C are 



(10) 



(11) 



(12) 



Cii ^ Ka2~ 3KRnlao/2 - Kklao + 771^ , 

C12 = -2KRnlT^/^kJj~k,nsTs, 

Ci3 = -2KRnsT!/^n'J-f~nX, 

C21 = kxUs , 

C22 = klusTs/j + "fns - 

C23 = k^TsTi'Jj ~ Kk^ai/A, 

C31 = -K > 

C32 = -Kk^ai/A, 

C33 = -{T^^'us + ryj/j + jUs-Kklao/i. (13) 



0-2 



i'q and 

The boundary conditions for the functions 0, 



(6) u and v are the following: 



9(0) = e'(l) = u'(0) = u'(l) = v{0) = v{l) = . (14) 

Equation (|l^) with the boundary conditions (^^ define a 
linear boundary-value problem: there are three boundary 
conditions at the base and three at the top. A simple nu- 
merical procedure (realized in MATLAB) enabled us to 
avoid the unpleasant shooting in three parameters. The 
procedure employs the linearity of the problem. We first 
complement the three boundary conditions at the base 
by three arbitrary boundary conditions at the base, and 
compute numerically three independent solutions of Eqs. 
( pO| ) . The general solution can be represented as a linear 
combination of these three independent solutions that in- 
cludes three arbitrary coefficients. Demanding that the 
three remaining boundary conditions at the top be satis- 
fied, we obtain three homogeneous linear algebraic equa- 
tions for the coefficients. A nontrivial solution requires 
that the determinant vanish, which yields the eigenvalue 
7. Varying R at fixed F, K and k^, we determine the 
critical value R = Rc for instability from the condition 
Re 7 = 0. We found this algorithm to be accurate and 
efficient. 

In the whole region of the parameter space that we 
explored we found that Im7 = at the instability on- 
set. Therefore, thermal granular convection does not 
exhibit overstability and can be analyzed in terms of 
marginal stability. Figure |^a shows the marginal sta- 
bility curves versus the horizontal wave number k^ at a 
fixed K and two different values of F The curves 

exhibit minima {k*,R*), similarly to the convection in 
classical fluids [Q. Therefore, the convection threshold 
in the horizontally infinite layer is R — R^. Close to the 
onset, the expected horizontal size of the convection cell 
is 2n/k*. Figure § depicts these convection cells obtained 
by plotting the field lines of the respective velocity field 
(u(jj) sin k^x , v{y) cos k^x) , found numerically. One can 
see that at small F the cells occupy the whole layer of 
granular gas and are elongated in the horizontal direc- 
tion. At large F the cells are effectively located near the 
base, and their aspect ratio is close to unity. 

Figure 2b corresponds to a zero gravity: _F = 0. Here a 
different symmetry-breaking instability occurs: the one 
that leads to phase separation |l^, [l^, |2l|, p2| . 
When R exceeds the marginal stability threshold R* {F — 
0), the laterally symmetric stripe of enhanced particle 
density at the top wall becomes unstable and gives way 
to a 2D steady state. In contrast to the convection, the 
new steady state with a broken translational symmetry is 
static: no mean flow. The quantity R*{F = 0) can be cal- 
culated analytically . In our present notation, it is de- 
termined from the algebraic equation coth /i = ^, where 
/i = (i?*/2)i/2. This yields r{{F = 0) = 2.8785 . . .. The 
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FIG. 2: The critical values of the relative heat loss parame- 
ter R for the convection instability (non-zero gravity, a) and 
phase-separation instability (zero gravity, b) versus the hori- 
zontal wave number k^. The Knudsen number K = 0.02. 




FIG. 3: The convection cells at the instability onset k = 
and R — R^. The Knudsen number is K — 0.02. Figure 
a corresponds to the solid curve of Fig. Here F = 0.1, 
k* = 1.8 and R'^ = 0.49. Figure b corresponds to the large-_F 
limit, when the granulate is localized at the base. Here F = 5, 
k% = 19, and _R* = 2.86. 

minimum of the marginal stability curve occurs here at 
/cj = 0, that is for an infinitely long wavelength. 

We found that the crossover between the two instabil- 
ities is continuous. The dependences of i?* and fc* on 
the Froude number F are shown in Fig. ^. One can see 
that the Rl{F) dependence is non-monotonic. A stronger 
gravity is favorable for convection at very small F (as 
Fig. ||a also shows). However, this tendency is reversed 
at ~ 0.16, and R* starts to grow with F until it satu- 
rates at large F. In its turn, k* goes down monotonically 
with F and vanishes at F = 0. The decrease is quite slow 
at intermediate F, but becomes very rapid at very small 
F. As the phase separation instability does not exist in 
classical fluid, this low-i^ behavior is unique for granular 
fluid. 

The large-F limit deserves a special attention. Here 
the granulate is localized at the base. This regime is 
convenient in experiment, as particle collisions with the 
top wall (which are in reality inelastic) are avoided. A 
natural unit of distance in this regime is A = To/mg, 

1 /2 

while the time should be scaled to A/Tq . Correspond- 
ingly, (n) is deflned now as the total number of particles 
per unit length in the horizontal direction, divided by A. 
After rescaling the Froude number F drops out of Eqs. 
(0)"® enters the problem only via the top wall po- 
sition H/X = F. For F ^ 1 the top wall can be safely 



moved to inflnity, and F drops out of the problem com- 
pletely. Therefore, at large F, the convection threshold 
R* should depend only on K. Our numerical results fully 
support this prediction, see Fig. ^. How should k* be- 
have at large F7 Let us reintroduce for a moment the 
"physical" (dimensional) horizontal wave number kxph- 
The scaled critical wave number fc* = k^^j^H = k*^i^\F. 
As the product k*pf^X is the scaled wave number in the 
newly rescaled variables, it should be independent of F 
at large F. Therefore, k* should be proportional to F. 
Figure ^ shows that the quantity k*/F indeed approaches 
a constant (that depends on K) at large F. 
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FIG. 4: The convection threshold versus the Froude num- 
ber F. The solid line shows the Rc{F) curve at K = 0.02. 
At large F the threshold _R* approaches 2.85.... The dashed 
line is the Schwarzschild's curve Rs{F) that gives a neces- 
sary condition for convection. The small-F asymptotics of 
the Schwarzschild's curve is Rs — F/2. In the limit of large 
F, Rs{F) approaches 1.06514.... The inset shows that, as 
F ^ 0, R* approaches 2.8785 . . ., the threshold of the phase- 
separation instability. 
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FIG. 5: The convection threshold versus the Froude number 
F. Shown is the k*{F) dependence at K = 0.02. The inset 
shows that k* goes to zero as F — » 0. For large F, k%/F — > 
3.80.... 

These results give, for flxed values of F and the 
necessary and sufficient criterion for convection. It is of- 
ten useful to also have a simpler and easier-to-interpret 
criterion, even if approximate. A simplified criterion for 
convection can be obtained by neglecting the viscosity 
and heat conduction terms in the linearized equations 
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®"®' t'^^t is by taking the hmit K ^ 0. As the vis- 
cosity and heat conduction act against convection, this 
procedure obviously yields a necessary, but not sufficient, 
criterion for convection. 

Without the dissipative terms, Eqs. coincide 
with the linearized equations of ideal hydrodynamics of 
classical fluid with specified static profiles of tempera- 
ture Ts{y) and density ns{y). Even for this idealized 
problem, the exact criterion for convection can be ob- 
tained only numerically, and the result depends explic- 
itly on the specific profiles Ts{y) and ns{y). There is, 
however, a simple and general limit here in terms of the 
Schwarzschild's criterion j^] that yields a lower bound 
for the convection threshold. The Schwarzschild's crite- 
rion guarantees that there is no convection if the entropy 
of the fluid in the static state S{ns,Ts) grows with the 
height, that is S'{ns,Ts) > for any y. For the nearly 
elastic hard sphere model in 2D the granular entropy in 
the dilute limit is S{n,T) = ln(T/n). [Importantly, we 
are not making here any additional assumption: this sim- 
ple constitutive relation for S{n,T) was already used in 
Eqs. Therefore, the Schwarzschild's criterion 

can be rewritten in terms of the static temperature and 
density profiles and their first derivatives. For a given F, 
the static profiles are determined solely by the relative 
heat loss parameter R. Therefore, the Schwarzschild's 
criterion yields a critical value R = Rs{F) so that at 
R < Rs{F) there is no convection. The opposite in- 
equality R > Rs (F) yields a necessary (but of course not 
sufficient) criterion for convection. How to find Rs{F)7 
At small enough R the spatial derivative of the entropy, 
S'{y), is positive at any height y > 0. By increasing R 
we observe that, at the critical value R — Rs{F), the 
entropy derivative S'{y) vanishes at some point y. It is 
crucial that this point is always y = 0. Increasing R fur- 
ther, we would already have an interval of heights where 
S'{y) < 0. Therefore, the Schwarzschild's curve Rs{F) 
can be obtained from the condition S'{y = 0) = 0. This 
curve, obtained numerically, is shown by the dashed line 
in Fig. IJ. As expected, the exact (necessary and suf- 
ficient) convection threshold curve always lies above the 
Schwarzschild's curve. 

The small-i^ and large-_F asymptotics of the 
Schwarzschild's curve can be obtained analytically. Let 
us first consider the case of F ^ 1. As will be seen 
from the result, here R 1 too, and one can repre- 
sent the steady-state solutions as Ts{y) = 1 + dTs{y), 
and risiy) = 1 + Sns{y), where STs ^ 1 and Srig -C 1. 
Substituting these expressions into Eq. (^) and keep- 
ing only the first order quantities, we obtain two very 
simple linear differential equations. Solving them with 
the respective boundary and normalization conditions, 
we obtain Ts{y) ~ I - Ry + (i?/2)j/^ and ns{y) ~ 
I + F/2 - i?/3 + {R- F)y - (i?/2)j/2. The condition 
S"(t/ = 0) ==0 then yields the desired small- asymp- 
totics: Rs{F < 1) ~ F/2. 



At F ^ 1, one can conveniently use the analytic solu- 
tion for the static profiles in the Lagrangian mass coor- 
dinate 1^. In this limit 



Ts{z) 



ns{z) 



liiVR/2)z 
/r72/2(z) ' 



(15) 



where z = (i?/2)^/^(l — n), and /„(. . .) is the modified 
Bessel function of the first kind. (Recall that here we 
rescale the Eulerian coordinate y by A = To/mg and 
define (n) as the number of particles per unit length in 
the horizontal direction divided by A.) Using Eqs. (|T^), 
we compute the y-derivative of the entropy: 



S' 



(16) 



This expression vanishes when Iq{z) — Azli{z) = which 
occurs at z = z* = 0.72977 .... As y = corresponds to 
/i = 0, we immediately obtain Rs — Izl — 1.06514 . . .. 
At R < Rs we have S" > everywhere, so R = Rs 
indeed corresponds to the Schwarzschild's criterion. The 
value of Rs{F ^ 1) = 1.06514 . . . shows up as the large- 
F plateau of the Schwarzschild's curve in Fig. ^ 

Now let us return to the exact (necessary and suffi- 
cient) criterion for convection, found by solving the full 
linearized problem numerically. The dependence of the 
convection threshold i?* on the Knudsen number K (at 
a fixed F) is shown in Fig. ^. R* grows with K, because 
the viscosity and heat conduction (both of which scale 
like K) tend to suppress convection ||]. As if the 
viscosity and heat conduction become negligible. Still, 
one should expect that the limiting value of the criti- 
cal heat loss parameter Rk^o {F) — linix^o ^* {F, K) 
is greater than the Schwarzschild's lower bound Rs{F). 
This is indeed what is seen in Fig. |6[ where the 
Schwarzschild's value Rs{F = 0.1) ~ 0.05 is shown by 
the empty circle. 




FIG. 6: Convection threshold R* versus the Knudsen number 
K a.t F = 0.1. As K ^ 0, the value of R* is greater than 
the Schwarzschild's value Rs{F = 0.1) ~ 0.05, shown by the 
empty circle. 
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SUMMARY 

We performed a linear stability analysis of the static 
state in a horizontal layer of granular gas driven from be- 
low. The hydrodynamic theory Q that wc employed in 
our analysis is expected to be valid when the mean free 
path of the particles is much smaller than any length scale 
described hydrodynamically. We have found the convec- 
tion threshold, in terms of the relative heat loss number 
_R, versus the two other scaled numbers of the problem: 
the Froude number F and the Knudsen number K. We 
have predicted the morphology of convection cells at the 
onset of convection. As F — + 0, the convection insta- 
bility goes over continuously into the phase-separation 
instability |l7[ |l|, |o[ |2^, |2|l . At large F the convec- 
tion threshold depends only on K. We established a sim- 
ple connection between thermal granular convection and 
classical thermal convection of ideal compressible fluid. 
The connection is given in terms of the Schwarzschild's 
criterion, a universal necessary (but not sufficient) con- 
dition for thermal convection. A further development of 
the theory should account for the excluded-volume ef- 
fects jl^. Importantly, the simple Schwarzschild's crite- 
rion will be readily available in the finite-density theory. 
Indeed, this criterion requires only the knowledge of the 
static profiles of the granular temperature and density 
and the constitutive relation for the granular entropy. 
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